{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Historical ice thickness\n",
    "\n",
    "This notebook plots ICESat-2 sea ice thickness vs. Pan-Arctic Ice Ocean Modeling and Assimilation System [(PIOMAS)](http://psc.apl.uw.edu/research/projects/arctic-sea-ice-volume-anomaly/) sea ice thickness. PIOMAS is a project out of the Polar Science Center at the University of Washington that models historical sea ice thickness and volume. Plotting ICESat-2 and PIOMAS sea ice thickness allows for a comparison of model predictions and satellite observations. \n",
    " \n",
    "**Input**:\n",
    " - xarray dataset from the jupyter book's google bucket\n",
    " \n",
    " \n",
    " **Output**: \n",
    "  - Plots of sea ice thickness "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{tip}\n",
    "Try running this notebook in Google Colab! Toggle over the rocketship icon at the top of the page and click Colab to open a new window and run the notebook. <br><br>To run a single cell, type **Shift+Enter**. To run the whole notebook, under **Runtime** click **Run all**. Note that you will have to run the notebook from the very beginning and load all the Google Colab dependencies for the code to work.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "#this cell will load dependencies for running the notebook in Google Colab\n",
    "#this cell may take a while to run\n",
    "import sys\n",
    "\n",
    "#if code is running in google colab, run these cells to install neccessary libraries\n",
    "if 'google.colab' in sys.modules: \n",
    "    !pip install netcdf4\n",
    "    !pip install xarray==0.16.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import notebook dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt\n",
    "from textwrap import wrap\n",
    "\n",
    "#remove warnings to improve display\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    " \n",
    "#increase resolution for notebook outputs\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data into notebook\n",
    "Copy file from the book's google bucket and load into an xarray dataset. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load ICESat-2 dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "!gsutil -m cp gs://is2-pso-seaice/icesat2-book-data.nc ./\n",
    "IS2dataset = xr.open_dataset('icesat2-book-data.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load PIOMAS dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-output"
    ]
   },
   "outputs": [],
   "source": [
    "!gsutil -m cp gs://is2-pso-seaice/piomas-regridded-data.nc ./\n",
    "PIOdataset = xr.open_dataset('piomas-regridded-data.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restrict datasets to region of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "#function from regional_analysis notebook\n",
    "def restrictRegionally(dataset, regionKeyList): \n",
    "    \"\"\"Restrict dataset to input regions.\n",
    "    \n",
    "    Args: \n",
    "        dataset (xr Dataset): dataset generated by Load_IS2 notebook\n",
    "        regionKeyList (list): list of region keys to restrict data to \n",
    "        \n",
    "    Returns: \n",
    "        regionalDataset (xr Dataset): dataset with restricted data to input regions\n",
    "    \"\"\"\n",
    "    \n",
    "    def checkKeys(regionKeyList, regionTbl): \n",
    "        \"\"\"Check that regionKeyList was defined correctly\n",
    "\n",
    "        Raises: \n",
    "            ValueError if regionKeyList was not defined correctly \n",
    "            warning if all data was removed from the dataset\n",
    "        \"\"\"\n",
    "        if type(regionKeyList) != list: #raise a ValueError if regionKeyList is not a list \n",
    "            raise ValueError('regionKeyList needs to be a list. \\nFor example, if you want to restrict data to the Beaufort Sea, define regionKeyList = [13]')\n",
    "\n",
    "        for key in regionKeyList: \n",
    "            if key not in list(regionTbl['key']): \n",
    "                raise ValueError('Region key ' + str(key) + ' does not exist in region mask. \\n Redefine regionKeyList with key numbers from table')\n",
    "\n",
    "        if len(regionKeyList) == 0: \n",
    "            warnings.warn('You removed all the data from the dataset. Are you sure you wanted to do this? \\n If not, make sure the list regionKeyList is not empty and try again. \\n If you intended to keep data from all regions, set regionKeyList = list(tbl[\\\"key\\\"])')\n",
    " \n",
    "    #create a table of keys and labels\n",
    "    regionMask = dataset.region_mask.attrs\n",
    "    regionTbl = pd.DataFrame({'key': regionMask['keys'], 'label': regionMask['labels']})\n",
    "    \n",
    "    #call function to check if regionKeyList was defined correctly\n",
    "    checkKeys(regionKeyList, regionTbl)\n",
    "    \n",
    "    #keys to remove (all keys that are note listed in regionKeyList)\n",
    "    keysToRemove = [key for key in list(regionTbl['key']) if key not in regionKeyList]\n",
    "    \n",
    "    #filter elements from the ice thickness DataArray where the region is the desired region\n",
    "    regionalDataset = dataset.copy()\n",
    "    for var in dataset.data_vars: \n",
    "        if var != 'seaice_conc_monthly_cdr':\n",
    "            regionalVar = regionalDataset[var]\n",
    "            for key in keysToRemove: \n",
    "                regionalVar = regionalVar.where(regionalVar['region_mask'] != key)\n",
    "            regionalDataset[var] = regionalVar\n",
    "    \n",
    "    #find name of labels \n",
    "    labels = [regionTbl[regionTbl['key'] == key]['label'].item() for key in regionKeyList]\n",
    "    \n",
    "    #add new attributes describing changes made to the dataset\n",
    "    if len(labels) < len(regionTbl['key']): \n",
    "        if set(regionKeyList) == set([10,11,12,13,15]): #convert to sets so unordered lists are compared\n",
    "            regionalDataset.attrs['regions with data'] = 'Inner Arctic'\n",
    "        else:    \n",
    "            regionalDataset.attrs['regions with data'] = ('%s' % ', '.join(map(str, labels)))\n",
    "        print('Regions selected: ' + regionalDataset.attrs['regions with data'])\n",
    "    else: \n",
    "        regionalDataset.attrs['regions with data'] = 'All'\n",
    "        print('Regions selected: All \\nNo regions will be removed')\n",
    "    \n",
    "    return regionalDataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myKeys = [10,11,12,13,15] #Inner Arctic\n",
    "IS2dataset = restrictRegionally(IS2dataset, regionKeyList = myKeys)\n",
    "PIOdataset = restrictRegionally(PIOdataset, regionKeyList = myKeys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate ICESat-2 means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate mean monthly ICESat-2 sea ice thickness \n",
    "IS2means = IS2dataset['ice_thickness_filled'].mean(dim = ['x','y'], skipna = True)\n",
    "IS2means.attrs['long_name'] = 'ICESat-2 mean ice thickness'\n",
    "\n",
    "#calculate mean monthly ICESat-2 sea ice thickness uncertainty\n",
    "IS2uncs = IS2dataset['ice_thickness_unc_filled'].mean(dim = ['x','y'], skipna = True)\n",
    "IS2uncs.attrs['long_name'] = 'ICESat-2 mean uncertainty'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compile into one dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IS2 = xr.Dataset(data_vars = {'mean_ice_thickness': IS2means, 'mean_ice_thickness_unc': IS2uncs})\n",
    "IS2.attrs = IS2dataset.attrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate PIOMAS means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean sea ice thickness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PIOmeans = PIOdataset['PIOMAS_ice_thickness'].mean(dim = ['x','y'], skipna = True)\n",
    "PIOmeans.attrs['long_name'] = 'PIOMAS mean ice thickness'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean detrended uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define functions used in calculation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_varDT(Years, Thickness):\n",
    "    \"\"\" Detrend linear time series  \n",
    "    \"\"\"\n",
    "    trendT, interceptT, r_valsT, probT, stderrT = scipy.stats.linregress(Years, Thickness)\n",
    "    lineT = (trendT * Years) + interceptT\n",
    "    ThicknessDT = Thickness - lineT\n",
    "    return ThicknessDT, lineT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_mon(xarray_val, desired_mon):\n",
    "    \"\"\" Check if current month is desired month \n",
    "    \"\"\"\n",
    "    xarray_mon = pd.to_datetime(xarray_val).month\n",
    "    if xarray_mon == desired_mon:\n",
    "        return xarray_val\n",
    "    else: \n",
    "        return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Call functions to perform uncertainty calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#loop through months in a year and calculate uncertainty for each month \n",
    "standDevDT, standDev, monDates, thicknessDT, slopes = [],[],[],[],[]\n",
    "months = np.arange(1, 12 + 1, 1)  #integers corresponding to months of the year\n",
    "\n",
    "for mon in months:\n",
    "    #get dates for each month (i.e. for January get 1978-01, 1979-01, 1980-01, etc)\n",
    "    times_by_mon = [is_mon(val, mon) for val in PIOmeans.time.values if (is_mon(val, mon) != None)]\n",
    "    monDates.append(times_by_mon)\n",
    "    \n",
    "    #get array with year-mon converted to integer value\n",
    "    time_integers = np.arange(0, len(times_by_mon))\n",
    "    \n",
    "    #get detrended uncertainty and append to list\n",
    "    thickness_detrended, regLine = get_varDT(time_integers, PIOmeans.sel(time = times_by_mon).values)\n",
    "    thicknessDT.append(thickness_detrended)\n",
    "    standDevDT.append(np.std(thickness_detrended))\n",
    "    slopes.append(regLine)\n",
    "    \n",
    "    #get trended uncertainty and append to list \n",
    "    standDev.append(np.std(PIOmeans.sel(time = times_by_mon).values))\n",
    "    \n",
    "#get a list as long as PIOmeans of all the monthly natural variabilities\n",
    "standDevCyclic = []\n",
    "for date in PIOmeans.time.values:\n",
    "    standDevCyclic.append(standDevDT[pd.to_datetime(date).month - 1])\n",
    "    \n",
    "#create DataArray with corresponding uncertainties\n",
    "PIOuncs = xr.DataArray(standDevCyclic, dims = 'time', coords = {'time': PIOmeans.time.values}, attrs = {'long_name':'PIOMAS mean detrended uncertainty'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compile into one dataset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PIO = xr.Dataset(data_vars = {'mean_ice_thickness': PIOmeans, 'mean_ice_thickness_unc': PIOuncs})\n",
    "PIO.attrs = PIOdataset.attrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Winter comparison: PIOMAS vs. ICESat-2\n",
    "Here we will plot PIOMAS and ICESat-2 mean monthly ice thickness on the same axes to visualize differences in the two datasets over winter 2018-2019 and winter 2019-2020. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define plotting function "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotMeans(PIO, IS2, title = None, figPath = None):\n",
    "    \"\"\"Plot ICESat-2 ice thickness and PIOMAS ice thickness on the same axis\n",
    "    \n",
    "    Args: \n",
    "        PIO (xr Dataset): dataset with PIOMAS mean ice thickness and mean ice thickness detrended uncertainty\n",
    "        IS2 (xr Dataset): dataset with ICESat-2 mean ice thickness and mean ice thickness uncertainty\n",
    "        title (str, optional): string description of plot (default to None)\n",
    "        figPath (str, optional): path to save fig (default to None)\n",
    "\n",
    "    Returns: \n",
    "        Figure displayed in notebook \n",
    "        \n",
    "    Notes: \n",
    "        Uncertainties commented out of function\n",
    "    \"\"\"\n",
    "    \n",
    "    #initialize figure and axes\n",
    "    fig = plt.figure(figsize = [8,5]) \n",
    "    ax = plt.axes()\n",
    "    gridlines = plt.grid(b = True, linestyle = '--', alpha = 0.4) #add gridlines \n",
    "\n",
    "    #plot PIOMAS data \n",
    "    PIO['mean_ice_thickness'].plot(ax = ax, color = 'purple', marker = 'o', label = PIO['mean_ice_thickness'].attrs['long_name'])\n",
    "    #ax.fill_between(PIO.time.values, PIO['mean_ice_thickness'] + PIO['mean_ice_thickness_unc'], PIO['mean_ice_thickness'] - PIO['mean_ice_thickness_unc'], facecolor = 'purple', alpha = 0.1, edgecolor = 'none', label = PIO['mean_ice_thickness_unc'].attrs['long_name'])\n",
    "\n",
    "    #plot ICESat-2 data\n",
    "    IS2['mean_ice_thickness'].plot(color = 'dodgerblue', marker = 's', linestyle = '-.', label = IS2['mean_ice_thickness'].attrs['long_name'])\n",
    "    #ax.fill_between(IS2.time.values, IS2['mean_ice_thickness'] + IS2['mean_ice_thickness_unc'], IS2['mean_ice_thickness'] - IS2['mean_ice_thickness_unc'], facecolor = 'dodgerblue', alpha = 0.1, edgecolor = 'none', label = IS2['mean_ice_thickness_unc'].attrs['long_name'])\n",
    "\n",
    "    #display plot with descriptive legend, ylabel, and title\n",
    "    plt.legend()\n",
    "    plt.ylabel('Mean sea ice thickness (m)')\n",
    "    plt.title(PIO.attrs['regions with data'] + ': mean monthly sea ice thickness')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot comparison for winter 18-19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMeans(PIO = PIO.sel(time = slice('Nov 2018','Apr 2019')), IS2 = IS2.sel(time = slice('Nov 2018','Apr 2019')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot comparison for winter 19-20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMeans(PIO = PIO.sel(time = slice('Nov 2019','Apr 2020')), IS2 = IS2.sel(time = slice('Nov 2019','Apr 2020')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot PIOMAS with overlayed ICESat-2 data\n",
    "From the PIOMAS reanalysis data, this plot shows the yearly melting and growing of sea ice as well as the long term decreasing trend in mean ice thickness. ICESat-2 data is overlayed on the PIOMAS plot, along with a rolling 12 month mean. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize figure and axes \n",
    "fig = plt.figure(figsize = [8,5]) \n",
    "ax = plt.axes()\n",
    "gridlines = plt.grid(b = True, linestyle = '--', alpha = 0.4) #add gridlines \n",
    "\n",
    "#plot data \n",
    "PIO['mean_ice_thickness'].plot(ax = ax, color = 'steelblue', label = PIO['mean_ice_thickness'].attrs['long_name'])\n",
    "IS2['mean_ice_thickness'].plot(ax = ax, color = 'magenta', label = IS2['mean_ice_thickness'].attrs['long_name'])\n",
    "\n",
    "#using pandas rolling mean function, add a rolling mean \n",
    "rollingMeanInterval = 12 #number of months over which to take mean\n",
    "tbl = pd.DataFrame({'means': PIO['mean_ice_thickness'].values})\n",
    "rollingMean = np.array(tbl.rolling(rollingMeanInterval).mean()['means'])\n",
    "ax.plot(PIO.time.values, rollingMean, c = 'darkblue', label = 'PIOMAS ' + str(rollingMeanInterval) + ' month mean')\n",
    "\n",
    "#display plot with descriptive legend, ylabel, and title\n",
    "plt.legend()\n",
    "plt.ylabel('Mean sea ice thickness (m)')\n",
    "plt.title(PIO.attrs['regions with data'] + ': mean monthly sea ice thickness')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Understanding detrended uncertainty \n",
    "Detrended uncertainty captures natural variability in mean ice thickness, while ignoring the downward trend of ice thickness due to climate change. The regression line is subtracted from the mean ice thickness resulting in detrended ice thickness, which can be used to solve for natural variability (detrended uncertainty) in ice thickness.<br><br>The figure below shows these values for a single month over the period of 1978-2020 to illustrate this concept. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choose month of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "month = 7 #July"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize trended & detrended thickness "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize figure & axes\n",
    "fig = plt.figure(figsize = [8,5]) \n",
    "ax = plt.axes()\n",
    "gridlines = plt.grid(b = True, linestyle = '--', alpha = 0.4) #add gridlines \n",
    "years = np.arange(1978, 1978 + len(thicknessDT[month - 1]), 1)\n",
    "plt.title(PIO.attrs['regions with data'] + ': ' + pd.to_datetime(monDates[month - 1])[0].strftime('%B') + ' ice thickness & uncertainty', fontsize = 14)\n",
    "\n",
    "#plot thickness \n",
    "monthlyThickness = PIO.sel(time = monDates[month - 1]).mean_ice_thickness.values\n",
    "ax.plot(years, monthlyThickness, color = 'blue', label = 'mean ice thickness')\n",
    "\n",
    "#plot regression line\n",
    "ax.plot(years, slopes[month - 1], color = 'magenta', linestyle = '-.', label = 'regression line')\n",
    "\n",
    "#plot detrended thickness \n",
    "ax.plot(years, thicknessDT[month - 1], color = 'maroon', label = 'mean detrended ice thickness')\n",
    "\n",
    "#add text describing variability \n",
    "explanation = \"\\n\".join(wrap('Observing the difference in spread of the two plots, one can note that natural variability from the detrended mean sea ice thickness is ' + str(round(standDevDT[month - 1],2)) + ' meters, while mean monthly variability (including the decreasing trend in ice thickness due to climate change) is ' + str(round(standDev[month - 1],2)) + ' meters.', 40))\n",
    "ax.text(x = 1.075, y = 0.73, s = explanation, horizontalalignment = 'left', verticalalignment = 'top', fontsize = 11.5, transform = ax.transAxes, bbox = {'facecolor':'none', 'edgecolor':'lightgrey', 'boxstyle':'round'})\n",
    "\n",
    "#add legends & show plot\n",
    "plt.legend(bbox_to_anchor=(1.56, 1), fontsize = 11.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
